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ABST RACT 

The 'Isotherm Wigration Method’ has teeu fou.nd to he 
an effoctive ni?,merical technique for the solution of the heat 
conduction problems wi bh phase cha,nga involving temperature 
dependent properties. The principal equation for this method 
is derived with the iso'thorm position, x(0,t) in the conventio: 
heat conduction equation. This interchange of variables 
minimizes the computer time quite a bit comparod to the other 
numerical methods. This is because the properties are evaluat 
only once at each isotherm temperature and are used throughout 
the computation with no need to recalculate them. In the 
present study, the divided difference form of the governing 
equation has boon modified to incorporate the unequal temporat 
intervals in ’ tiravo-tomperature’ grids so that the exact proper 
ties ob Gained oxperimontaliy at unequal temperature intervals 
can bo used in the computation without interpolation. Moreove 
usually the phase change temperature does not fall on one of 
the regular isotherms taken for convenience at round temperatu 
of 10°, 20°, 30° etc. or 100°, 200°, 300° etc. Henco this 
modified IMM enables one to use any random temperature values. 

In the present work, the modified IMM has been described 
to determine the temperature profiles and the interface 
location in terms of X(i,j) for unidimensional solidification 
of tho saturated and unsaturated liquid under various geometri 
and different boundary conditions involving temperature 



dcpendont properties too. Tho insido aiad outside freezings 
for cylindrical/splierical geometries have also been calculated 
by this method. The system chosen to study the method is 
pure water at 0°C (for the saturated case) and at 20'^C (for 
the unsaturated case) . The validity of the technique has 
boon checked with good accuracy in certain cases available 
in the literature. 



[^1 


iTOMI'lITClATURE 
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© . K 
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— for B.C. of the 2nd kind in cylindrical/spherical 


^^c ’ ^ geome try 



9 - & 


for B.C. of the 3rd Icind in the solid zone 


op 

00 

9-0 

= Q — for B.C.s of the 1st and 3rd kind in the 
“ ^0 liquid zone 

T(i) Dimensionless temperature of the ith isotherm 
TIQ; Defined in Bq.(32) 

t time 

X dimensionless position of the isotherm 
X(i,j) the position of the ith isotherm at jth time 
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CHAPTER 1 


IHTRODUCTIOEf 

1. 1 The Importance of Solidificai;! on /Melting Problems ; 

The change of state occurring with melting or 
freezing is associated with many of today's industrially 
relevant problems. The heat transfer problan with phase 
change is evident in the solidification of casting, freezing, 
thawing of soils and foodstuffs, preservation of medical 
and biochemical materials, thermal storage devices, the 
ablation of missile skins under aerodynamics, heating, 
welding etc. Because of the economic importance of such 
processes, a wide range of research in the effects of phase 
change and in the physical models for it has been conducted. 

!• 2 Difficulties in Establishing an Exact Solution for the 
Solidification /Belting Problems ; 

In the transient heat conduction problem, the 
supply or removal of heat from any element causes it to 
change its temperature. In the case of solidification or 
melting, however, the system is essentially comnlicated 
due to the release or absorption of the latent heat of 
fusion which brings the non-linearity in the governing 


diff erential equations. 
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1 • 5 Different Approaches to Solve tlie Phase Change Problems ; 

A variety of approaches to the solution of phase 
change problems in solidification or melting classified as 
’Stefan’s Problems' have been used by different investigators. 
Solutions have been obtained for various geometries and boundary 
conditions for well over hundred years. After the reviews by 
Bankoff [12] and Muehlbauer et al. [13], over 400 papers have 
appeared in the International Journals [14]. Since exact 
analysis is difficult for all but very simple problans, 
various approximations have been used by different investigators. 
These namely are conformal mapping, embedding technique, 
finite elanent method, heat balance integral tec'i-inique, 
integral equation method, iterative analytical method, n\mierical 
methods, perturbation method etc. Most of the above methods 
assume constant thermophysical properties for the system 
which is far from the reality. In some cases linear or 
other kinds of variations in properties with temperature are 
assumed. These are generally solved by the conventional 
nomerical method, where space and time grids are used. The 
properties are evaluated at each grid noint due to the 
change in temperature at the grids with the passage of time. 

This results in excess computer time usage. 

1, 4 A New Numerical Approach ’Isotherm Migration Method’ (IMM) ; 

An effective way out has been' suggested [l, 2] to 
transform the conventional heat conduction equation so that 
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its solutions express space as a function of the temperature 
and time coordinates. Instead of writing temperature 
© = 0-(x,t), one expresses 'x' as a function of temperature 
and time i.c. x=x(0-,t). Then, using conventional numerical 
technique the positions of the various isotherms at different 
times are calculated. Hence one essentially traces the move- 
ment of the isotherms (i.e, constant tanperature lines) at 
different times. So it is called the ISOTHEHVI MIGRATION 
METHOD (IMM), This obviates the need to calculate the 
properties at each grid point. The values of the properties 
once calculated, can be used throughout the computation with 
no need to recalculate them. This process minimizes the 
computer time quite a bit. The efficacy of this technique 
over other approximate methods lies in the heat transfer and 
diffusion problems with property variation. Apart from 
this, in ablation analysis the effect of surface recession 
is more easily treated. The temperature distribution is more 
readily applicable to thermal stress analysis since more points 
occur in the region of large temperature gradients usually 
at the specimen surfaces, 

1.5 Literature Review on IMM ; 

Ohemousko [l] has explored this idea for the 
solution of non-linear heat conduction problems in a medium 
with phase change, Dix and Gizek [2] used this technique 
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for transient heat conduction analvsn s in one dimension for 
the boundary condition of the t^'ird kind. Crank and Phahle [3] 
have applied this method to solve the problem of the melting 
of a plane sheet of ice. Crank and Cupta [5] have extended 
the IMM to two dimension system evaluating the problem posed 
by the solidification of a square prism of fluid initially 
at constant temperature. 

1.6 The Modification of the Eiciating IMjVI in the Present Study 
and the Advantage Iherein ; 

In the present study the IMM has been modified to 
include unequal temperature grids in the finite difference 
form for the numerical calculations. This modification 
eliminates certain shortcomings in the existing IMM, The 
existing formulation of the IM. requires the use of constant 
taittporaaure intervals which often lead to interpolation between 
values of properties obtained experimentally at specific 
temperatures. More serious is the fact that the phase change 
tonperature does not usually fall on one of the regular 
isothems taken conveniently at roimd temperature of 10°, 20°, 

30° etc. or 100°, 200°, 500° etc. Hence, this modification 
helps to use the exact properties at unequal temperature 
intervals reported in the literature, without interpolation. 
This also cuts down the ccmputer time a bit and gives less 
erroneous results. This method can co'nveniently be extended 
to solve multiphase solidification problems where interface 
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temperatures usually are not at regiilar intervals. 

1. 7 Pocus on the Details of the Present Work ; 

It is known that a solution to a -Droblein in solidi- 
fication is also a solution to the corresponding inverse 
prohlons in melting, the present work has been arbi-'^rarily 
restricted to the process of solidification in order to 
facilitate the detailed search. 

The modified IMM has been used to calculate the 
interface location and the temperature profiles as a function 
of time for a undimensional solidifying system subjected to 
any possible combination of the following conditions: 

1, Solidification can take place in Cartesian, 
cylindrical or spherical geometries. 

2, The boundary condition at the surface of the 
body can be of tho 1st, 2nd or the 3rd kind. 

3, In the case of cylindrical or spherical body 
the solidification can proceed outward or 
inward. 

4, Initial state of the liquid could be at the 

freezing temperature 0-^ (saturated case) or 

at any temperature 0- above the fusion temperature 

s 

(unsaturated case). 

5, The property variation as a function of temperature 
in the system (both in the solid and liquid 
region) can also be taken into account. 
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6, The temperature intervals can be equally or 
xmoqually spaced in the timo-t emp erature grid. 

The system treated in each caso is further defined by the 
following postulates; (a) The phase change occurs at a fixed 
temperature, the fusion tannerature (©-q). (t) Convection 
effects on heat transfer due to densitv chapge has been 
ignored, (c) Initially the liquid is at a uniform temperature, 
(d) In the constant property and unsaturated case, the 
property of solid and liq-uid has separately been evaluated 
at the fusion temperature • 



GHAPTER 2 


MATHMATICAL F0MULATI0R3 

2,1 The Transfonned Ea’uation for IMM and the Significance of 
Each Term Involved ; 

The heat conduction equation in one dimension system 
without internal heat generation is expressed as 






( 1 ) 


^.n S X 

Here n is called the shape factor, 
n = 0 planar geometry 
n = 1 cylindrical gcomet2y 
n = 2 spherical geometry 

Bq. (l) can be transformed to the desired form to 
suit IMM with the help of following relations; 

( 2 ) 


(3) 


(— ) - 

( Bxs-l 


and 



' -a t ^x 

= - 

/ X 

t'' 2>t 

Substituting Eqs 

( 2 ) and 

(3) in 

O n i ( 1^9 ) ( ^ ] 

^ 1 



1 


“ 1 

Fe 


If 




< ^X n-1, 


.n' 


S /''B X- 

■■BO 


-1 




8 


'M 


-II 






f/AlL .-ri,_T7- „ ..n-l ^'^Xw x\-2 


ji 




x""+K.n.x" 






) 


1 

Oarcelling (>~3’)"’ from both sides of the equation, we ha-we 


or 


or 


- (|f’ 


(|t) 




^E j , K n r ^2e^V B X\-l x ^ -2 


■c)Q 




Xx ( 


) 


K 


( 


/^C 


w vX \-2 n 1 /D K \ / \-l I 

>1 W q; - X “ K J 


(|f) = X = a 


W AEs. ^ ^ X x -2 n -1 /'^x -1 1 

~ K ^ c)'9'''‘c>©^ 5 


(4) 


Thus Sq. (4) becomes the governing equation expressed in 
terms of x(©, t) for the one dimensional unsteady state heat 
conduction. The left hand side of eq, (4), indeed represents 
the migration rate of a particular isotherm (©) as a function 
of time. The first terai on the right hand side of Eq. (4) 
implies the migration rate of the isotherm through the 
geometry whose surface area is not a function of thickness 
(i.e. planox* system) and the therraonhysical properties are 
independent of temperature. The second term takes into 
consideration of the shaoe factor to include the cylindrical 
and the Spherical geometries. The last term accounts for 
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the temperature dependence oi a and K. The above Eq. (4) 
has already been derived from energy ba3.ance [2] which 
provides better insight into the physical situation and 
allows easier incorporation of the other eff.ects, for example, 
phase change and temperature dependent heat generation. 


2. 2 Mfferent Boundary Conditions Studied for the Solidification 
Problem : 

The IMM has been used for a \midimensional solidifying 
system whose surface is exposed to the following boundary 
conditions. 


(a) The boundary condition of the 1st kind is 


0 -r 0,^ at x-O 

W 


t 0 


(5) 


(b) Tho boundary condition of the 2nd kind is 

1o = ^ If = 

(c) The boundary condition of the 3rd kind is 


K 




= ) at x=0, ■*» ^0 


W CD 




( 6 ) 


(7) 


W 00' 

The condition at the solid-liquid interface describing tho 
process of solidification is 

„ _ p X(.^) 

^l^cje-'l ^2''&0^2 “ o ^ ^ t'' 

at x=x^ 


( 8 ) 
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The Gqiiation is valid ior the solidification of the 
•unsaturated liquid initially at the temperature (© ), 

For the saturated liquid Bq. (8) simplifies to 


tr -Pi 


at x=x 

o 

The initial conditions are defined as 

© = 0g for the unsaturated liquid 


and 


© = ©. 


for the saturated liquid 


(9) 


( 10 ) 

( 11 ) 


2.3 Non-dimensionalisa.tion of the Transfomed Equation and 
the Boundary Oonditions Involved ; 

To generalise the application of the IMM so that any 
phase change problem can easily be treated with the one 
governing equation, it is required to non-dimensionalise 
Eq. (4) and also the boundary conditions involved. The choice 
of variables to define the dimensionless parameter is made 
according to the boundary conditions and shapes describing 
the solidification syston. 


(a) Boundary Condition of the 1st kind, 

© - © 


© - © 

m S£ 

■^1 - © - © » 

0 -W 


rn .8 

2 - Ss- 


n 

= a t/L , for planar system 


= t/R^, for cylindrical/spherical system 
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Oy 


and X = x/l for nlanar s‘'''steiii 

= r/R for cylindrical/sphorical system 
Sub sti tilting thu above dimensionless parameter, the Eq. (4) 
takes the foim 


ki _ « ff C>^ X w'c)X'-2 

V't " %L 


n 

X 


1 

K 




?)X n-1 

~b 


J 

(12) 


The dimensionless Eq. (12) gives the location of different 
isothclms both in the solid and the liquid region as daioted 
by the subscripts 1 and 2 respectively. 

The Eq. (8), representing the interface location, takes the 
form 




[Zi(-^)l.(oyV = 


at X = (13) 

The Bq. (l3) is applicable for the unsaturated liquid ^ while 
for the saturated liquid the Eq. (9) takes the form 


D> 


(^) _ ^o^^o '^w^ / ^Xn~1 


(14) 


The boundary condition at X=0 and t ^ 0 


T = 0 


(15) 
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(Td) The bomdary condition of the 2nd kind. 

Here the dimensionless parameters are 

for plan'X system 

for cylindrical/sphorical system 

p 

i: = t/1 for planar system 

p 

= t/R for cylindrical/spherical system 

and X = x/L for planar system 

= r/R for cylindrical/spherical system 
With these parameters Eq. (4) transforms again into Eq. (l2), 
whereas Eqs (8) and (9) become 


T = 


K • Q 
o 


a. E 


and 




and 


X 


1o ■ ^ 


^ “o 

" fo*-% 


r^_ ^ V Vz^-: 


(. 

T^l K T 

0 w J- Q 


h -^1 




+ 


K 

•IT 


2 


] 


[At the interface = K^] 


(Ai.) 

''at-' 


a X a/ ^ 1 

0 o 


( 16 ) 


(17) 


at X=X^ • 
o . 

f or nnsaturated and saturated liquid respectively. The 
boundary condition at X=0, and i.e, at the surface 

represented by the Eq. (6), becomes 



(c) The "boundary condition of the 3rd kind. 


Here the dimensionless parameters chosen are: 



a 


®-co 


O’ - 0- ’ 

o oo 


T, 




t: = t/1 for planar system 

= oc t/H for cylindrical/spherical system 
X = x/1 for planc-r system 

= r/R for cylindrical/spherical system. 

Hero also Eq. (4) tran^orms into eq. (12)^ wnereas Iqs. (8) and 
(9) "become 



2>Xn~1 

£>T^1 








!3 Z ' 
■^1^2 


and : 

~ A 

at X = X^ 

The "boundary condition 
represented by Eq.(7), 


00^ ^o / ^X^-1 

Ta 

rospectivoly. 

at X=0 and x^O, 
takes the form. 


(19) 


( 20 ) 


i.e. at the surface 



h, Il» , m 
K 


( 21 ) 


The initial conditions which are same for all the B.C.s are 
T 2 = 1 for un saturated liquid; T^=l for saturated liquid (22) 
The numerical solution of Bq#(l2) with different boundary 
conditions has been dealt with in Chapter 3» 



lUMaRICAL OALCULATIOHS 


At first, it is observed that the transfoimed Bq. (4) 
or Bq. (12) is non-linear even for the constant thermal 
properties and Cartesian system, whereas the original heat 
conduction equation is linear. This is relatively ■unimportant 
if an explicit difference method is used as in Bq.(25), but 
would call for the solution of a set of non-linear algebraic 
eq'uations if an implicit difference method is introduced. 

Hence explicit or forward difference method [10] is used here. 
The non-linearity causes no problon in applying this method. 

3,1 Divided Difference Form for the Governing Equation and 
the Boundary Conditions ; 

To incorporate the unequal temperature interval in 
(t-T) grid for the numerical solution, the finite difference 
form for the 2nd order derivative is expressed in the 
following form: 



If the temperature interval is equally spaced, i.e, T(i+l)-T(i) 

= T(i)-T(i-l) = j^lT, then Eq, (23) simplifies to the conventional 
form 

® j Z(i+l,j)- 2X(i,j) + X(i-l, 3 ) 

■ T^ (A T)^ 


( 24 ) 
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So, forward difference in time and central difference in 
temperature yield an equation from, tne Bq. (l2) for the ith 
isotherm position at the time (3+I) as follows: 




The Bq.(25) can be siiUTvlified to the form 

- C, n C, HTI 

X(i,j+l) 5!^A+XB) •C^+(l-C2) •X(i,3 ) - 

( 26 ) 

where 

XA = X{i+l,3)*(T(i)-T(i-l)) (27) 

XB = X(i-1,3)- (1(1+1) *l(i)) (28) 


Q ^ 0: (1 ) 
1 ~ 

o 


, g ( l .( m .L " . l (,i.rl)..) ) 

(X(i-k-l, 3 ) “X (i-l, 3 ) ) ( T( i+1 )-T(i )> ( T(i )-T ( i-1 ) 


(29) 

(30) 


Cg = 0 ^ (T(i-k‘l)-T(i-l)) 
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1 cdll 


ft. 


At^ 


TICI 


K(i+1) - Z(i-l) 


(31) 

(32) 


XI = (X(i+l,j) - X(i~l,j)) 


(33) 


Similar^ the finite difference form for the botmdary conditions 
at the interface and- at the wall for different cases becomes. 

(a) For B.C. of the 1st kind 

At I = X 

o 

xd.d^l) = X(i,3) ^ 

(34) 

For saturated liquid it becomes 

.w ' , . (& - (4 ) • (T(i)-T(i-l) ) 

-(i.3+1) =X(i,j) - [ (X°i,3V A(i-i,3)) ] (35) 

At X-0 , t ^ 0 

Xj(l, 3 ),= 0 for planer system 

X(l,j) =1,0 for cylindrical/spherical system 

(b) For B.C. of the 2nd kind 

At X=Zq The Eqs(l6) and (l7) become 


A 'C'q.-l 

X(i,3+1) = X(i,3) + 


T(i)-T(i-l) ^ 


-1 


^X(i+l,j) - X(i,j)^_3^_^ 
- k 7 ^(i+1) - T(i) ^ ^ 


( 36 ) 
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and 

X(i, o+l) 


Up I<)-(T( i )-T(i -1)) 
(^o ^ )-(X(i,j)-X(i-l,r)') 


(57) 


At the surface i.e. X=:0 and t 0 , the corresponding 
Bq. (l8) becomes 

\=[t(2) -X(2,o+ 1)J /K(i)A;^ (38) 


Eq.(58) is for planer system, wiiile for cylindrical and 
spherical system the Eq. (l8) takes the form; 

\ = ri(2) + (p*(l-X(2,j+l)))}A^(i)Ao (59) 

Here p=l for inside solidification and p=-l for the outside 
solidificati on. 


(c) For B.C. of the 5rd kind 

At X=X^, Bqs.(l6) and (17) for the unsaturated and saturated 
liquid become 


X(i,D+l) = X(i, j) + 






-1 


-i 




(40) 


and 

X(i,3+1) = X(i,j) 


A-c 0^i%- &-^)*(T(i)-T(i-l)) 
(X(i, j ) - X(i-1, j ) ) 

(41) 


At the surface of the body i.e, X=0, ^ 0, the finite 

difference forn of Eq. (19) becomes. 
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= 1(?)/(T + for rlaner s'^T-gtem 

(42) 

= 1 ( 2 ) / ( 1+p • ( " ^ ' £ 't~l ) ) f 0 r cylindri cal/ 

spherical system 

(43) 

Here also p=l for inside freezing and p= -1 for outside 
freezing. 

3.2 Convergence. Stability Criterion and Truncation Error ; 

The explicit schane necessitates the convergence and 
the stability criterion to be satisfied to impls-ment the 
numerical solution. The convergence implies that the finite 
difference approximation will approach the exact solution 
when the size of the increments employed is made infinitesimally 
small. The stability implies that errors associated with 
the use of increments of finite size or round off error will 
not grow as the calculation proceeds. When finite difference 
procedure is both convergent and stable, a comparison of the 
calculations made using two different increment intervals is 
generally a good indication of the reliability that raa'.?- be 
assigned to the results. Although some definite tests are 
available for the convergence and stability criterion cf the 
linear equations, but, no straight forward treatment to the 
convergence and the stability of the non-linear equations is 
readily available. In this situation, one must proceed with 
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the physical concept of the formulations. 

The stability criterion for non-linear equations 
can be approximated in the following manner; 

(1) The non-linear equation, if possible, is simplified 
to the linear form, just by eliminating some parameters or 
variables involved. Then the stability criterion for the 

simplified linear equation is found, low small variation 
or perturbation is intuitively given to the deteimined value 
to take into accoimt the non-linear eff ect/eff ects of the 
main equation, finally one should accept that perturbated 
value of the stability criterion for which convergence 
criterion is achieved, 

(2) If it is not possible to simplify into the 
linear form, then, intuitively satisfying the physical concept 
of the problem, one can establish a relation from the coeffi- 
cients of the divided difference fom of the non-linear 
equation. 

The g-oveming Eq. (l2) belongs to the second category, 
Eq.(l2), even for the simplest case like planer geometry 
and constant property, is non-linear and its divided difference 
form i.e. Eq. (26) simplifies to the expression 

X(i,j+l) = (XA + XB)-C^+(l- (44) 

Uow for X(i,3+l) to increase with an increase in 
X(i,j), as expected physically the coefficient of X(i,3) must 
be positive i,. e., 
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(1 - C 2 ) 0 


O'. [1 


a(l) 
a 


( 




o (Ki+l, j)-X(i~l,3))2(T(i+l)-.T(i))*'(T(i)-T(i-l)) 


)] 


This finally yields 


> 0 


< 




( 45 ) 


S5.nce the equation is non-linear a more exact analysis is 
not possible. For the equal temperature interval, the 
expression (45) simplifies to the form as given in [2] 

AX < (X(i-H,;j) - X(i-l,o))^/8.(a(i)/aQ) (46) 

or in the dimensional form 


At < (x(i+l,3) - x(i-l, j ) )^/8, a (i) (47) 


From a comparison of tho Taylor scries representation of 
Eq. ( 12 ) and its differential equivalent, the truncation error 
is fo\md to be nroportional to (at) and (AT) which compare 
with and ( X)^ for the es^nlicit equation. 

3.5 Comnutational Procedure for Different Oases ; 

It is evident from Eq, ( 26 ) and the boundary conditions 
at X=0 and X=X , tihat if X(i,j) is known at one time t, one 
can compute the X(i,o) at (x^ Ax) time i.e. X(i,j+l). This 
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corresponds "bo all isotherms chosen including the interface one. 
This implies that the present state of the system can be evaluated 
by this explicit technique if the past state of the system is 
known. Thus one must generate an initial set of isotherms to 
begin tho If®!, Analytic or some alternative solution is needed 
to provide the temperature distribution at some small time from 
which IMM can proceed on. In the present study the analytical 
solutions [4,9 and 15] have been used, according to the different 
cases to initialise the isotherm locations including the interface 
isotherm for small time t=l min. In certain cases this IMM has 
been used to compute the solidification even upto 10 hours. The 
results thus obtained are found to be in good agreement with [4]. 
For clarity in the graphs only tho computations upto 1 hour have 
been plotted in all tho cases. 

For very short times the ci^iindrical and spherical solutions 
approach tho planar solution. Accordingly, the results of the 
planar solution [4] for t=l min have been used to start the IMM 
in the cylindrical and spherical cases. The small error from this 
approximation dies out rapidly with time as was confirmed by the 
calculation done with different initial temperature distribution. 

For the boundary condition of the 2nd and 3rd kind, the 
rate of surface temperature increase is determined from an energy 
balance near the surface as represented in Eqs. (38,39,42 and 43). 


The first appearaiice' of the new isothorm has been located inside 
the body near th# wall by linear approximation with those of 
other existing isotherms. It does no-^ affect the resul-^ because 




A 
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tlivj tumpui’aturj intorvalG choson aro small enough like -5°, -4°, 
-3^,-2*^, -I® ,0®C so on. A.s the su.rfaco temperature decroasos, 
tho- new isotherms must be created at the boundary and allowed to 
migrate into the modium. 

The time increment has been calculated satisfying the 
stability criterion from the Eq.(45). To calculate from Eq. 
(45), the two alternate isotherms nearer to the surface havo 
been taken because they are closely spaced in comparison to other 
distant isotherms. After every five iterations js-rr has been 
recalculated from the Eq.(45) since the spatial distance b>itweun 
any two isotherms increases as the timo progresses. This pro- 
cedure also reduces the amount of computation time. The effects 
of 3rd and 4th terms in Bq.(26) for csrlindrical/spherical geometry 
and property variation respectively have been taken into account, 

choosing conveniently ^x = 0,'15l A rl (^j^2.culated* 
because their contribution in bh'- main Eq.(26) is small for the 
system studied. It has been found that in case of inside soli- 
dification, for cylinder or sphere, the computation becomes 
unstable even before reaching the neighbourhood of the centre. 
Since X(i,j) becomes lower, the '^d term in Eq. (26) dominates 
and coasos to give the right solution. To combat this situation 
has been cut short accordingly in different stages of 
computations so that solidification can be achieved upto the 
neighbourhood of the centre. The complete solidication hrus been 
calculated by extrapolating the situations of partial solidi- 
fications upto the neighbourhood of the centre. 



CHAPTER 4 


PLOVf DIAGRM OP THE COMPUTER PROGRAMMIHG 
The outline of the computer programming comprising of 
all the system variables studied in the present work has been 
shown below in the form of block diagrams. The listing of 
the final FORTRAN IV program is given in the Appendix II. 

Codes used in the program and their significances 
with appropriate values; 

(a) Geometry types (N) 

= o Planar geometry 
= 1 Cylindrical geometry 
= 2 Spherical geometry 

(b) Boundary conditions (ROD) 

= -1 Boundary condition of the 1st kind 

= 0 Boundary condition of the 2nd kind 

= 1 Boundary condition of the 3rd kind 

(c) Initial state of the liquid (NIS) 

= 1 Solidification of saturated liquid 

= 2 Solidification of unsaturated liquid 

(d) Temperature dependent property (VIP) 

= 1 Constant properigr 

= 2 Variable properly 

(e) Inside/Outside Solidification (INOU) 

= -1 Outside freezing in cylindrical/spherical 

= 1 Inside freezing in cylindrical/spherical 

geometry 
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50,/'^ 

- " qiOD ;■ 


CalculatG the inter- 
face loca.tion from 
Sqs(54 or 35) 


I Calculate the inte 
'face location from 
|Sq.(36or 37) and 
the surface temp's- 
lerature from 
Eg. 


Calculate the inter- 
face location from 
Eg. (40 or 41) and the 
surface temperature 
from Eg. (43) 


Calculate the intei 
face location from 
the Eg. (34 or 35) 


After every 

five 

iterations , 

calculat 

t:ie At from eg. 

(45) 



Calculate the inter- 
face location from 
Eg. (36 or 37) and 
the surface tempera' 
ture from 7:lg.(38) 


Calculate the inter- 
face location from Eg, 
(40 or 4l) and the 
surface temperature 
from Eg. (42) 


=t -f t 


V ^ 


>.( 40 ) 


isotherm forme? 


1- Yes 


E+1 

lW+1 


(STO^ 


- -y . 

Calculate the position of the newly 
formed isotherm by the linear approxi- 
mation of the previous two isotherms 
and also its physical properties 


lOO 








GjiAP'JER S 


RESULTS MD DiaOTTflOTnT? 

Computations following the modified IMII with suitahl^r 
Ciiosen At from LQ..(45^) were performed for the selidi^icat; on 
of water initially at 0°C T saturated case] and at 20°C 
[unsaturated case] on DJC 1090 computer. The boundary- 
conditions chosen arbi-brarily, are the followings s 

(1) 0^ = -50°C for B.C. of the 1st kind. 

(2) qH = 0.1 cals/sec cm^ for E.C. of the 2nd kind 

■^x=o 

(3) h = 100 Btu/hr.ft2°f for B.C. of 'the 3rd kind 

= 4.18598 X 10~^ca,ls/sec.cm^ °C 

The physical propeities of water and ice used in the study have 
been reported in the Appendix I. Results were obtained in the 
form of A(i,;j) as a function of time. They have been shown 
graphically for certain cases and compared with the available 
literature values. The legend for each of the cases has been 
labelled in the figure itself, 

5.1 Results of the IMM for Solidification of Saturated liquid 
Under Various Conditions ; 

Pigs. (1 to 3) compare the IMM with the analytical 
solutions [4] showing the interface location and the temperature 
profiles for the Cartesion solidifying systems subjected to 
the boundary cnnditions of the 1st, 2nd and 3rd kind respectively. 
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It ia io'infi liha.''; for the J.C, of the 1st kind, the isotherms 
are closel;;!,' sn-'rced at small times and the spatial distance 
between any two isotherms increases with the passage of time. 
For the B.C. of the 2nd id.nd (Fig. 2), the isotherms are more 
or less equally spaced thj-oughout the time of solidification, 
since constant amount of heat is drawn all along the process 
Jrom the surface; while for the B.C. of the 3rd kind (Fig. 3) 
the spatial distance between any two consecutive isotherms 
increases slowly with time. For a given flux (Fig. 2) and 
convection (Fig. 3) at the boundary, the surface temperature 
decreases with time which causes the inclusion of the new 
isotherms at the boundary which afterwards keep on migrating 
through the medium. 

Fig. 4 represents tlie comparison between the results 
of the IMl-1 for equal and unequal temperature interval in 
(T— t) grid. Here the temperature intervals are spaced at 
0°C, ~15°C, -25°0 and -30^0 for the numerical calculation 

done by the IMM. The response is found to be good. 

Fig. 5 depicts the interface location and the temperature 
profiles both for constant and the variable properties. In 
the constant property case, the properties are evaluated 
at the interface temperature. In the variable property case 
the initialisations of the IMM have been done from the modified 
error function solution [9]. The response follows the usual 
trend af the graph. Since the thermal conductivity of ice 
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1 iicroaGow vj.l bh decreasin^^ temperatiure [H] which facilitates 
the heat transfer through the iDody, both interface and 
other isotherms mi rate at a slightly faster rate than in 
the constant property case. 

Pig. 6 shows the comparison of the IMM with other 
approximato solutions [8] and [6] for depth of outward solidi- 
fication through cylindrical and spherical bodies respectively. 
The results agree well. Here solidification is faster in 
case of the cylinder than in the case of the sphere. This 
is so because the increment of surface area resulting from 
solidification is more rapid in case of the cylinder than the 
sphere. This process facilitates the higher heat transfer 
and thereby higher migration rate of the interface in case 
of the cylinder. 

For the inward solidification as shown in Pig, 7, the 
effect of solidification on the surface area acts in the 
reverse way, wherein the sphere freezes faster than the 
cylinder. Here the complete solidification has been shown. 
This has been achieved by the extrapolation of the partial 
solidification calculated by the IMM upto the neighbourhood 
of the centre. 

5.2 Be suit 3 of the IMM for Solidification of Unsaturated 
Liquid Under Various Conditions ; 

In Fig, 8 the res\xlts of the IMM have been compared 
with the Neumann's solution [4], Here temperature profiles 
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have been calculated both in the solid and the liquid regions. 

The location of the last isotherm corresponding to 20°C in 
the liquid region is calculated by the extrapolation technique^ 
Sjince no definite relation or boundary condition can be 
assigned to that point. Ideally, the temperature in the 
infinite medium asymptotically approaches its initial value 
T as 'X' becomes arbitrarily large. However, the concept 
of 'penetration distance' is used to locate the last isotherm 
in the liquid medium, Hor this purpose the following artifice 
is suggested. A parabola is fitted through X(i,j) of 10° 0 
isotherm and X(i,j) of 15° C isotherm ’with dT/dX = 0 at T=Tg 
ie at 20°C isotherm. The distance at the point of tangency 
is taken as the penetration distance. Its value is given by 

X(i~l,j) - vlT'(i-ZT^T(Trr7(f a-l)-f(i-2)) 

( 1 , j ) = ^ ) Jt (i ) ) / ( T ( 

Since the concept of the 'penetration distance' is not so 
well defined, the results of the IMM deviates from the analytical 
solutions [4] within 5 per cent for distant isotherms in the 
liquid region. However, in the solid zone the values are in 
good agreement with the Heumann's solution [4]. 

In Fig. (9), the results of the IMM for constant and 
variable properties with boundary condition of the 2nd kind 
have been plotted. Temperature fields are drawn both in the 
solid and in the liquid zone. The nature of the plots is 
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same as that of l''ig*(2). In tho unsaturated case, the 
interface aiid other isotoorus migrate through the medium at 
a slower rate incompari&on to the saturated case. The ItM 
hq.(26) has been initialised from Boley’s solution [151 for 
small time t=l min. Since Boley's solution is based on 
finite geometry with one side insulated, if cannot be used to 
compare the results of the lim for higher time. 

The temperature profiles and the interface location 
for temperature dependent property and unequal temperature 
interval in a planer solidifying system have been shown in 
tho Pig. (lO) . Comparison has also been made wit]i equal 
temperature interval case. The resTilts, thus obtained, are 
reasonably closed to each other. 

Pig. (11) depicts tho results of the IMM computed by the 
equal and the unequal temperature grids in tho divided 
difference form. In unequal temperature interval case, the 
isotherms are randomly chosen at -20° C, -15°C, -10° C, 0°C, 
-1°C, -3°C, -5°C, -8°C etc. Here the system chosen is 
cylindrical body exposed to the boundary condition of the 
3rd kind. The nature of the graph is satisfactory and the 
results in both the cases are reasonably close', to each other. 

Fig, (12) shows the effect of property variation with 
temperature occurring in the spherical geometry. The tempera- 
ture profiles in the solid and liquid zones and the position 
have been compared with the constant property case with the 
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propertiec evaluated at the interface temperature. The 
gradual increase of the thermal conductivity in the variable 
property case along the direction of heat flow facilitates 
the heat transfer and consequently the isotherms move at some- 


what faster rate. 
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Fig.l: Comparison of the results of the IMM with the 

analytical solutions [4] for boundary condition 
of the 1st kind 
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OOHC LUSIOl 

The above results in most of the cases compare with 
the literature values within 2 per cent. Tho deviation is 
slightly more (within 5 per cent) in the liquid region when 
the unsaturated liquid solidifies. To improve upon this 
one should look for better extrapolation techniques for 
locating the initial temperature isotherm. Convection in 
the liquid is ignored. However 5 , the major effects of variation 
of the thermal conductivities and heat capacities have been 
studied. Hence, the analysis presented here permits a ! 

realistic evaluation of tho interface location and temperature I 
fields for unidimonsion solidiCication incorporating both 
the equal and unequal temperature interval in tho time- I 

1 

temperature grids. | 

Though it is not a very major advance, it does increase [ 

( 

the utility of the IHM quite a bit. Concentration dependent , 

I 

diffusivity problems in mass transfer operations can be ' 

« I 

solved similarly. Future work can be directed towards extension 
of the method to two and three dimension and tho development 

of implicit method of solution of the equations. This modified ’I 

1 

IMM will be one important numerical tool to facilitate the 
multiphase solidification problems involved in the casting of 
steel and other alloys. 



Based upon the discussion during the thesis defence, 
the following has been added. 

(a) Eadiation Boundary Condition: 

The boundary condition for radiation at the surface 
X=0 and t 0, can be expressed in the dimensionless form as 


where h^ is radiation heat transfer coefficient 

6 F. P (t"^) 

A e w^ 


o p(\ 4 


6 is Stefan-Boltzman constant, cals/sec cra^(^C) 
is view factor, tiere = 1 

Pg is emissivity factor = -j-- - j— 


Substituting all these in Eq..(A), we have 
^ i/b'x = — Sg. — 2 


Taking divided difference form, we have 

■ <5 Pp li A 

= 1(2) t2,J+l) 1" 


6 P^ I X(2,J+1) 


T‘^ + T - T(2) = 0 
w w ^ ^ 


T_. can be calculated from the expression (B) by Regular PcQ ci 

W 

Method, Once the surface temperature is known, the remaining 
procedure is same as that of other boundary conditions alrc.-iy 
discussed. 
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APaJI^fDlX I 

Tho Physical Properties of Ic c j and Water and Their Temperature 

Dopo ndences 

(a) The required thermophysical properties of ice and water 
at 0°C [11] are listed below: 


Properties 


Water 

Units 

Donsib3)-( -f ) 

0.92 

0.99986 

gms/ cm^ 

Specific heat (C) 

0.492 

1.0092 

cals/gm °C 

Thermal conductivity (K) 

49.9205 

13.25x10“^ 

cals/sec cm 

xlO"^ 


°C 

Thermal diffusivity (a) 

0.0114 

1.321x10“^ 

cm^/sec 


The latent heat of solidification of water is 
^ = 79.78 cals/gm. 

(B) Thu variations of the specific heat of ice with temperature 
[11] are as follows; 

Tempera t\ire , Specific heat of ice, cals/gm °C 


0 

0.492 

-20 

0.4633 

-40 

0.4346 


Por small temperature range, it has been evaluated from the 
relation 

C = a + bt + ct^ where t is in °C 

a = 0.492; b = -5.5875x10”^ 
c = 1.793 X 10“^ 


and 
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(c) Tho variations of tlio specific heat of water with 
temperature [11] are as follows; 

Temperature, °C Specific heat, cals/gm °C 


0 

5 

10 

15 

20 


1.0092 

1.0051 

1.0021 

1.0001 

0.9983 


(d) The thermal conductivity of ice varies with temperatures [11] 
according to the relation 

I 

IC^ = Kq [1 - 17 X 10“^t] -170°C <" t < 0°C : 

whore t is in °C 

and K = 49.9205 x 10"'^ cals/sec cm.°C ; 

0 I 

(e) The variation of the thermal conductivity of water with 

temperatures has been reported in literature [11], J 

according to the relation 

^ ^20 [1+0.00281 (t-20)], 0°C <t < 80°C S 

I 

o i 

where t is in 0 i 

—*5 o 

and K 2 Q = 1.420 x 10 cals/sec. cm 0. 

IIo definite relation or discreate values for thermal diffuaivity ' 

1 

with temperature have been found in the literatures. So it , 

has been calculated from the relation a = K/ C, where K and C 
are evaluated at the required temperature from the above 


correlations. 
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COMPUTEH PROGRAMMING FOB 


isotherm^, MIGRATION .METHOD , CI,M,M) 


ONE DIMENSIONAL PHASE CHANGE PROBLKM(SOLIDIFICATION OF WATER} IN 
CARTESIAN.CTLINDRICAL AND SPHKRICAIr STSTEM, CONSTANT/ VARIABLE 
PROPERTIES.EQUAL/UNEQAL TEMPARATUKE INTERVALS.LIQlflD IS AT 
SATURATED/ONSATURATED STATE.THE SYSTEM IS SUBJECTED TO THE 
BOUNDARY CONDITIONS OF THE iST/2NU/3RD KIND, THE VALUE REQUIRED 
TO START THE ALGORITHIM IS FED BY HAND COMPUTATION FROM NEUMANS 
SOLUTIONS 




DEFINATION OF DIFFERENT PARAMETERS USED IN THE PROGRAM 




I 

J 

N 


INDEX 


K 


NIS 


TEMPERATURE GRID 
TIME GRID INDEX 
A FACTOR TO IDENTIFY THE COORDINATE SYSTEM,N=0 PLANER, 
N=1 CYLINDRICAL AND N=2 SPHERICAL CORDINATE 
TOTAL NUMBER OF ISOTHERM IN THE SOLID AND LIQUID 
ZONE INCLUCING THE INTERFACE ONE 
A FACTOR TO REPRESENT THE INITIAL STATE OF THE 


VIP 


LIQUID? NISal, SATURATED ?NIS=2, UNSATURATED LIQUID 
= A FACTOR TO INCORPORATE PROPERTY VARIATION WI 


:th 


INOU 


TEMP, V IP=1, CONST. PROP, ?VIP=2jr VARIABLE PORP, 

A FACTOR TO REPRESENT INWARO/OUTWARD SOLIDIFICATION 


KOD 


FOR CYLINDRICAL/SPHERICAL GEOMETRIES? INOU=l, INSIDE 
OR PLANER S0LIDIFIC.?IN0U=-1,0UTSIDE SOLIDIFIC, 

A FACTOR TO REPRESENT DIFFERENT BOUNDARY CONDITIONS 
KOD=-lilST, KIND?KOD«0, 2ND, KIND?KODBlt3RD, KIND 
A COUNTER TO SPECIFY THE INTERFACE ISOTHERM 
THERMAL DIFFUSIVITY OF THE SOLID AT THE 
INTERFACE TEMPARATURE 
ALPLO » THERMAL DIFFUSIVITY OF THE LIQUID AT THE 
INTERFACE TEMPARATURE 

ALPU) s THERMAL DIFFUSIVITY OF THE ITH ISOTHERM 


INF 

ALPSO 


RALPCDs ALPCIVALPSO 


TKSO 


OR TKH(I)/TKSO 


TKLQ 


TKHCI) 

RHO 

X(I,J) 

t(i5 


SPHS 

SPHL 

AL 

TAUMAX 

STIME 


DELTAU = 
CAT 


THERMAL CONDUCTIVITY OF THE SOLID AT THE 
interface TEMPERATURE 

THERMAL CONDUCTIVITY OF THE LIQUID AT THE 
INTERFACE TEMPERATURE 

THERMAL CONDUCTIVITY OF THE ITH ISOTHERM 

DENSITY OF SOLID OR LIQUID (ASSUMING CONSTANT) 

DIMENSIONLESS POSITION OF THE ITH ISOTHERM AT JTH TIME 

DIMENSIONLESS TEMPARATURE OF THE ITH ISOTHERM 

DEFINED IN THE TEXT 

SPECIFIC HEAT OF THE SOLID 

SPECIFIC HEAT OF THE LIQUID 

SUTABLY CHOSEN A PARAMETER IN DIMENSION OF LENGTH 
MAXIMUM TIME UP TO WHICH THE RESULT IS REQUIRED 
THE SMALL STARTING TIME AT WHICH THE ALGORITHM 
HAS BEEN INITIALISED 


H 

FLUX 



= HEAT TRANSFER COEFFICIENT OF THE SURROUNDINGS 


= CONSTANT HEAT FLUX DRAWN FROM^THE SURFACE 


If i(t ^ it; f f » :|c t i|| 


f 

r 

i 


OlOO 

C 

0200 

C 

)300 

C 

)400 

& 

1700 

C 

1800 

1900 

5 

000 

!ioo 

151 

i200 

1300 

152 

1400 

1500 

1600 

1700 

161 

1800 

|900 

10 

JOOO 

hoo 

111 

2200 

2300 

121 

2400 

2500 

2600 

131 

2700 

2800 

2900 

2000 

JlOO 

J200 

141 

3300 

1400 

20 

B500 

3600 

3700 

3800 

30 

3900 

40 

4000 

4100 

i200 

4300 

4400 

35 

4500 

4600 

42 

4700 

4800 

59 

4900 

5000 

5100 

5200 

5300 

43 

5400 

5500 

11 

5600 

100 

5700 

5800 

22 

5900 

200 

6000 

33 

6100 

6200 

6300 

300 


1000), TKHC45),S(4S),X(45, 1000), TKG(45) 

DfRAUFi<(&D} f G ALP 1 25 ) 

READCl ,5)K,NIS,INF,N,INOU,KOD 
FORMATCbIS) 

Tif PEI 5 1 

format (//,1 OX, 'VALUE OF DIFFERENT FACTORS USED IN THE PROGRAM') 

4< X Sfc 3f* M 

FORMAT C 1 0 X , ' *'*********"“*******"*"^»"»^»w«**^»«*w»w^»««*^«i«i*mwwf*(*»*i**»«*i*i*,«pi* ^ j 

TJfPE161,K,NlS,INF,N,INOU,KOD 

FORMAT C //, 2X,*K=', 13, 5X,*NISs', 13, 5X, ' INF*' , 13 , 5X, *Ns ' , 13 , 
25X,'IM0U«'I3,5X,*K0D=‘I3) ,aj, 

REApCl ,10)ALPSO,ALPLO,RHO,DELTAU,SPHS,SPHL,ALATEN,TKSO, 

5TKL0, VIP (CAT 
F0RMAT(bF9,6) 

TVPElll 

FORMATC//,10X, 'PHYSICAL PROPERTIES OF WATER AND ICE AT OC*) 

T I PL I i 1 

TJPE141 ,4l.ATEii,iKs6,TKLo7vfp?jAf^*' ' 

5F9"6’'^x'‘vit3‘'FPH4'F^®}?''' 

TIMjxi366o,0:iTiMEJl{oy6;Ai=i5,0 

?&y!!!$PS^‘"??9*l?J^^^''^^^**2);sfAUB(ALPS0*STIME)/CAL**2) 

IF vNlD«Ly» 1 ) tv-sii^F 

READ(l,20lCT(I),I=l,K) 

F0RMATC7F7,5) 

REAOCl ,30) U(I,i),I=l,K) 

F0RMAT(6F9,7) 

READCl,40)CTKH(I),Isl,K) 

IF(VIP,EQ,1)GDT035 
READ(l,lO)CRALPCl),Isl,K) 

F0RMATC6F10,8) 

H=4*18595E-03;TINF=0,0;FLUX=0,l 
IF CK0D,EQ,1)DTEM=0, 03333 
IF (ROD. EQ.O)DTEMaO, 00333 


F9.6,2X,'RH0a»,F9,6,2X, 


IFCVIP 


- .*,ltO)RALPCl)=l,0 
TlPE42e(TCl) ,l»lfK) 

F0RMAT(//,2X, ’SPECIFIED ISOTHERMS ARE VlX , 1 1 (2X, F7 .5 ) ) 
TyPE59 
FORMATCIX, 


TYPE43,(Xa,l),I»l,K) 
F0RMATC//,2X, 'INITIAL 
TYPESg 

IF(N,EQ,2)GOTOIOO 

if(n,eo.i5goto2oo 

TXPEll 

F0RMATC///,15X, 'PLANNER 
GOT0300 
yu pg2 2 

FaRMAT(///,15X, 'SPHERICAL 


POSITION OF THE ISOTHERMS' /IX, 1 1 ( IX, F9, 7) ) 


SOLUTIONS'///) 
SOLUTIONS'///) 


GOTD300 
TyPE33 

FORMATC///, 15X, 'CYLINDRICAL SOLUTIONS'///) 

J=1 ; JOal ; TKSOaO ,499205E«02 

KOUNTbT 

TAU»STAU+DELTAU?TL=STAU?TAUIM»0,180;NAT»0 





0100 
0200 
0300 
0400 
500 
600 
700 
800 
900 
000 
100 
200 
300 
400 
500 
600 
700 
800 
900 
000 
100 
200 
i|oo 
iioo 
ioo 
loo 
ioo 
00 
00 
•00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
§0 
00 
00 
00 
00 
00 
00 
0 
0 

00 
00 
00 
00 
00 
0 
0 

oo 

Is 
00 

DO 
0 
0 
|o 

p 

|o 
oo 
io 
m 
^0 
loo 


4000 

720 

1000 


CALCULATION OF THK ISOTHERM POSITION IN THE FROZEN AND 


D0400I=2,K 

IF(l.EQ,K)GOTO90O 

IFCK0D.GE,0)GQTU79 

XCl »J+i)sl,0 

IFCN.EO,0)XCl,J+l)=OtO 

IJ=J+i 

XI*(X(I + l,J)*XCI-ljj,J))^ ^ 

XA=XCI+1 , J)*{TCI)-TCI-1) ) 

XB=XCl-l.J)f(T(I+l)-T(I)) 

TI=CT(l+l)*T(I-l)) 

TTaCT(I+l)-T(l))+(TCI)-T(I-l)) 

TTI = C2,0=»*TI)/(TT*(XI**2) ) 

INF2=INF+1 

IF(i,EO.INF2)TKHCl-l)aTKLO 

TK1=CIKH(H-i5“TKHCI-1)) 

IFCVIP.EQ.ilfKIaO.O 
IFCVIP,EO.i •)RALPCI)»lfO 
IFCNIS.EQ.DGOTOSOO 
IFCl.EQ,INF)GOT0600 ^ 

IFa,Gf,INF,AND,VIP.EQ.l)RAl,P(I)=0,115877l 

TlsRALP?I)»6ELTilu»TTI 

T2ari*TI 

T3=RALP<I)*DELTAU 

XCIf J+1)»(T1 *TxA+X&) )+((!. *T 2)*XCI,J)>-C (T3*N)/X( I, J)) 
5-C(T3*TKl)/CTKHCI)»XI)) 

GOTU400 

SHL=(XCI+1#J)-X(I,J)5/(TCI+1)-TCI)) 

IF(K0D)76j.77i76 

calculation of the interface position FDR UNSATURATED CASE 


mmmmmmmmmmmm* 


XCI/ J+13=XCI »J)t(DELTAU)fC(TKSO*30,/SHS) + (TKLO=t‘20,/ 

5SHL) ) / ( ALPS04RH0*ALATEN) 

GOT0400 

XCI/J+l)»X(l£j)+(DELTAU#FLUX»AL)*((l,/SHS)-(TKL0/5HL*TKS0))/ 

SCALPSO^RHO^ALATEN) 

GOTO400 

CALCULilT?5N^OF^THE*^POSITlDN OF THE FARTHEST ISOTHERM IN THE 
UNFROZEN REGION BY PARABOLA FITTING METHOD 


TSKaSQRT(CTCl-l)»TCl))/CT(I-2)-TCI)n,_ ^ 

X(Ij-J + l)aCX(I-lrJ+l)«T5K*X(l»2,J+l))/(l«»T5K) 

G0T0400 

sH2=(xci,j)-x(i«i,a))/CTa)-T(i-i)> , 

CALCULATION OF THE INTERFACE POSITION FOR SATURATED CASE 
IF(X0D)82, 83,82 

X(I,J + l)aX(l,J) + CDELTAU>«'SPHS*30.0)/CSH2»ALATEN) 

G0TU400 

X ( I^ J + l )sX( I , J) + C FLUX*DELTAU^AL) / C ALPS0»RH0*ALATEN*SH2 ) 
IFCK0D)4000,14,13 

IFCJ,GE,7)GOT0710 , , ... 


X C I, U f § f / VIW X U / * V 

i?: '0E.Tw= 

iri CI&6-TC)tGE!!J?0l2]G0TO7 30)GOT0796 
TYPE1000,J,TAU,DELTAU, (X( I , J+ 1 ) , 1=1 ,K);TLsTAU 
GOT0796 


DELTAU»',F10,8,/1X,U 


820 

2000 

796 


AFTER EVERY ITERATION THE INCLUSION OF THE NEW ISOTHERM 
IS CHECKED FOR THE BOUNDARY CONDITIONS OF 2NDI.3RD KIND 


if(n;ge 

IFCN.EQ 

G0T086 

IFCN.GE 

IFCN.EO 


,l)TWN=CTC2)+INau»(l,-X(2,J+l)))/RAbP(l) 
,0)TWN=(TU) + INOU*(0,0-X(2,0 + 1)))/(TKH(1 


(TKHCn/TKSO) 


4.r v.’»«5i*UF’ACT=(Cl.-XC2.J+l}D»H»AL)/TKH(l) 
IF(N,EO,0)FACT=(X(2,J+U»H»AL)/TKH(l) 
TWN=r(2j/(l.+lN0U*FACT) 

IF(TWN,LE,TINF)GOT0202 

TDIFsTCU-TWN 

IF(N,E0.0)XCl,J+l)=XC2,J+i)=t‘((T(l)-TWN)/CT(2)-TWN)) 

IF(N.GE,l)Xa£J+l)=X,+(XC2,J+l)-l,)*((TCl)-TWN)/CT(2)-TWN)) 

DDIFk(DTEM-TDIF);IFC0DIF)39,39,49 

IF(J,GE,2)G0T0796 

G0TD820 

D041 1I=1,K 

SClDsTCIl) 

D051 IIal,K 
YCIi,J+UsX(II,J+U 
D0711I®! ,K 
GALPCII)=RALPCH) 

TKGtlDeTKHCII) 

KsK + 1 
INF-INF+1 

NEW ISOTHERM IS ADDED IN THE SYSTEM 


SYSTEM 


D0444I*2,K 

TCI)sS(I-n 

CONTINUE 

TCl)=TCl)»DTEM 

D0555Is2,K 

XCI#J+UaY(I«l,J+l) 

CONTINUE 

LOCATION, OF THE NEW 


ISOTHERM 


LINEAR APPROX, WITHIN 


SMALL TEMPERATURE INTERVAL 


IF ( N EQ 0 

IF(N!GEll)XCi;j+D»!l, + a(2,J+i)-l.)^C(TU5-f 
D0666l*ZtK ^ 

RALPCl)»GALP<I**n 

TKHCI)»TKGCI-1) 

CAL CULATIOM OF THE THERMAL CONDUCTIVITY AND 
DIFFUSIVITY OF THE NEW ISOTHERM ACCORDING TO 
APPROXIMATION AS REPORTED IN THE LITERATURE 


xa,j+i)ax(2{j+n*( 

XC1,J+i5»!1, + CX(2,J+ 


V-VMT, 


-TWN)) 


THERMAL 

LINEAR 


IF(VIP,GT, 


,KOD,|:Q,UTKHU)=TKSO*a,-0.g51*(T< 


IFCVIP,GT,l.0.AND,KOD.EQ,0)TKH(l)=TKS0*<l.»0,5l08*Ti 

IFCVIP,GT.l,)RALPU)=RALP(2)+(RALP(K-5)-RALPC2))*CT( 

5(T(K-5)-TC2)3 

if(vip,eq,i,o5tkhcu=tkso 

IF(VIP.EQ,1,03RALP(1)=1,0 

TYPE2000,KOUNT,a,TAU,DELTAU,TWN,(X(I,J), 1=1 ,K-n 

TYPE2000,KOUNT, IJ,TAU,DELTAU.TWN, (XCliJ+l ) , I«l - 

FORMAT(2X,I5,2X,I7,2X,Fi2,6,2X,FiO,8,2X,F9,5/ix,9(l 

JsJ+l ; TAU=TAU+DELTAU 

IFCTAU,GE,TAUMAX)GOT0202 

IFU,Ge,998)GOT0204 

aj=u«Jo 

IFCJJ,GE,10)GOT0203 

G0TQ302 


U-U)) 

i)iTC2))/ 


j,i=i ,R-i; 

X,F9,5/ix.9(lX,E13,8)) 


